Code
library(zoib)
library(tidyverse)
library(GGally)
library(ggpubr)
library(kableExtra)
library(mice)
library(parallel)
library(ggplot2)
library(glue)
theme_set(theme_pubr(legend = "bottom"))Objective 1
1 Methods
The primary question of interest is to determine if lower pre- or post-op attention scores are associated with a lower percentage of oral feeds at discharge. Because the outcome of percent oral feeds at discharge is bounded on the interval [0,1], a zero and one inflated beta regression model was a natural choice to model the relationship between attention scores and the outcome. Because the data only has 132 neonates and 24 potential covariates, a simple linear regression of the outcome on each potential covariate was fit to determine which covariates had significant linear influences on the outcome. This led to the choice to control for the length of time the neonate was intubated as well as its cardiac anatomy as covariates. Although not showing a significant linear influence on the outcome, age and gender of the neonate were included in the model as they are important demographic factors that are often controlled for in the literature. Lastly, we chose to run two separate models, one with pre-operation attention scores and the other with post-operation attention scores as covariates to answer the research question.
The zero-inflated model that was fit can be written as: \[ f(y_i|\eta_i) = \begin{cases} p_i, & y_i =0 \\ (1-p_i)q_i, & y_i =1 \\ (1 − p_i)(1 − q_i)Beta(\alpha_{i_1}, \alpha_{i_2}) &y_i \in (0, 1) \end{cases} \]
\[ \begin{aligned} logit\left(\mu^{(0,1)} = \frac{\alpha_1}{\alpha_1+ \alpha_2}\right) = \mathbf x \boldsymbol \beta_1 \\ log\left(v_i = \alpha_1+ \alpha_2\right) = \mathbf x \boldsymbol \beta_2 \\ logit\left(p_i\right) = \mathbf x \boldsymbol \beta_3 \\ logit\left(q_i\right) = \mathbf x \boldsymbol \beta_4 \end{aligned} \]
Where our response, \(Y_i\) represents the percentage of oral feeds at discharge of the \(i^{th}\) neonate, \(f\left(y_i|\eta_i\right)\) is the pdf of \(Y_i\), \(p_i\) is the probability of the \(i^{th}\) neonate having no oral feeds at discharge, \(q_i\) is the probability of the \(i^{th}\) neonate having 100% oral feeds at discharge, given \(y_i\neq 0\). \(Beta\left(\alpha_{i_1},\alpha_{i_2}\right)\) is the probability distribution of \(Y_i\) given \(y_i\neq 0\) and \(y_i \neq 1\). \(\mathbf{x}\) is a design matrix containing covariates representing the intercept term, attention score, length of intubation, cardiac anatomy, age, and gender.
The model was fit in a Bayesian framework using the zoib package in R. Diffuse normal priors, with a precision term of \(10^{-3}\) were set for all parameters. For each imputed data set, the models were run with four chains that yielded 1400 draws of the posterior distribution, after tuning parameters were set. We then pooled the chains together and checked autocorrelation and trace plots to ensure proper mixing and convergence. We then performed a check of the linearity assumption by plotting residuals against their fitted values. Lastly a posterior predictive distribution check was made by making nine replicates of the data and visually assessing differences between the distributions of the replicates and the data.
2 Results
When pooling the multiple imputation of the twenty datasets together, the models provide little to no evidence that a relationship exists between percent oral feed at discharge and the pre- or post-op attention scores. The results of both models, one including pre-op attention scores and the other including post-op attention scores, were exponentiated to be on the odds ratio scale and are shown in Table 1 and Table 2 respectively. For the dichotomous attention score models, there was evidence of a relationship between the the probability of having full oral feeds at discharge and pre-operation attention scores (see Table 3), yet there was no evidence of the association existing with post-operation scores (Table 4).
To check the assumption of linearity we plotted the residuals versus predicted outcome, which showed strong linearity (see Section 3.1.3). For a posterior predictive check, nine replicated samples were drawn from the posterior distribution and compared with the original data on percent oral feed at discharge (Figures Figure 1 and Figure 2, Figure 3, Figure 4). The histograms of the replicated samples do not reflect the clear zero-inflation of the original data. Further analysis in a thousand draws of the posterior distribution showed that the maximum or minimum (zero and one) were never drawn. A similar check conditioning on the outcome being between zero and one also yielded little confidence in conditional model fit. The poor model fit to our data suggests our results to be highly suspect. Future work should address these concerns, potentially by changing the covariates in the original model- such as adding a design matrix for modeling the log variance- or selecting another model entirely.
Code
nnns_imputed<- readRDS("../Haojia_work/nnns_imputed.rds")
dat <- lapply(1:20, function(i) complete(nnns_imputed, i))Code
# ZOIB MODELS
set.seed(11282023)
tictoc::tic()
# Define a function for model fitting
fit_model <- function(d) {
zoib(
Percent_of_feeds_taken_by_mouth_at_discharge ~
Pre_Op_NNNS_attention_score +
Length_of_intubation_days +
Cardiac_Anatomy +
Age_at_Surgery_days +
Female #x1 design matrix
| 1 | #x2 design matrix
Pre_Op_NNNS_attention_score +
Length_of_intubation_days +
Cardiac_Anatomy +
Age_at_Surgery_days +
Female | #X3 design matrix
Pre_Op_NNNS_attention_score +
Length_of_intubation_days +
Cardiac_Anatomy +
Age_at_Surgery_days +
Female, #x4 design matrix
data = d,
n.response = 1,
zero.inflation = TRUE,
one.inflation = TRUE,
link.mu = "logit",
link.x0 = "logit",
link.x1 = "logit",
random = 0,
n.chain = 4,
n.iter = 3000,
n.thin = 2,
n.burn = 200,
seeds = c(11, 29, 20, 23)
)
}
model_results <- lapply(dat, fit_model)
# Save results
saveRDS(model_results, "pre_op_models.rds")
tictoc::toc()
#Post-op Model
set.seed(11282023)
tictoc::tic()
# Define a function for model fitting
fit_model <- function(d) {
zoib(
Percent_of_feeds_taken_by_mouth_at_discharge ~
Post_Op_NNNS_attention_score +
Length_of_intubation_days +
Cardiac_Anatomy +
Age_at_Surgery_days +
Female #x1 design matrix
| 1 | #x2 design matrix
Post_Op_NNNS_attention_score +
Length_of_intubation_days +
Cardiac_Anatomy +
Age_at_Surgery_days +
Female | #X3 design matrix
Post_Op_NNNS_attention_score +
Length_of_intubation_days +
Cardiac_Anatomy +
Age_at_Surgery_days +
Female, #x4 design matrix
data = d,
n.response = 1,
zero.inflation = TRUE,
one.inflation = TRUE,
link.mu = "logit",
link.x0 = "logit",
link.x1 = "logit",
random = 0,
n.chain = 4,
n.iter = 3000,
n.thin = 2,
n.burn = 200,
seeds = c(11, 29, 20, 23)
)
}
# Parallelize model fitting
model_results <- lapply(dat, fit_model)
# Save results
saveRDS(model_results, "post_op_models.rds")
tictoc::toc()
post_op_models =readRDS(file="data/data/post_op_models.rds")
pre_op_models =readRDS(file="data/data/pre_op_models.rds")
post_op_coeff = list()
pre_op_coeff = list()
for(i in 1:length(post_op_models)){
pre_op_coeff[[i]] = pre_op_models[[i]]$coeff
post_op_coeff[[i]] = post_op_models[[i]]$coeff
}
pooled_pre_op = runjags::combine.mcmc(mcmc.objects = pre_op_coeff, collapse.chains = FALSE)
pooled_post_op = runjags::combine.mcmc(mcmc.objects = post_op_coeff, collapse.chains = FALSE)
saveRDS(pooled_pre_op, "pooled_pre_op.rds")
saveRDS(pooled_post_op, "pooled_post_op.rds")Code
pooled_pre_op =readRDS("data/pooled_pre_op.rds")
pooled_post_op=readRDS("data/pooled_post_op.rds")
summary_pooled_pre_op = pooled_pre_op |> summary()
summary_pooled_post_op = pooled_post_op |> summary()
rnames = c("Baseline","Attention Score","Length of Intubation (d)",
"Single Ventricle w/ Arch Obstruction","Two Ventricles w/ Arch Obstruction",
"Age", "Female")
pre_mean = summary_pooled_pre_op$statistics[,"Mean"]
pre_lb = summary_pooled_pre_op$quantiles[,"2.5%"]
pre_ub = summary_pooled_pre_op$quantiles[,"97.5%"]
post_mean = summary_pooled_post_op$statistics[,"Mean"]
post_lb = summary_pooled_post_op$quantiles[,"2.5%"]
post_ub = summary_pooled_post_op$quantiles[,"97.5%"]
pre_b_df= cbind("Mean"= pre_mean[1:7] |> exp() |> round(2),
"2.5%"= pre_lb[1:7] |> exp() |> round(2),
"97.5%" = pre_ub[1:7] |> exp() |> round(2))
pre_b0_df= cbind("Mean"= pre_mean[8:14] |> exp() |> round(2),
"2.5%"= pre_lb[8:14] |> exp() |> round(2),
"97.5%" = pre_ub[8:14] |> exp() |> round(2))
pre_b1_df= cbind("Mean"= pre_mean[15:21] |> exp() |> round(2),
"2.5%"= pre_lb[15:21] |> exp() |> round(2),
"97.5%" = pre_ub[15:21] |> exp() |> round(2))
pre_df = cbind(pre_b_df,pre_b0_df,pre_b1_df)
rownames(pre_df) = rnames
post_b_df= cbind("Mean"= post_mean[1:7] |> exp() |> round(2),
"2.5%"= post_lb[1:7] |> exp() |> round(2),
"97.5%" = post_ub[1:7] |> exp() |> round(2))
post_b0_df= cbind("Mean"= post_mean[8:14] |> exp() |> round(2),
"2.5%"= post_lb[8:14] |> exp() |> round(2),
"97.5%" = post_ub[8:14] |> exp() |> round(2))
post_b1_df= cbind("Mean "= post_mean[15:21] |> exp() |> round(2),
"2.5%"= post_lb[15:21] |> exp() |> round(2),
"97.5%" = post_ub[15:21] |> exp() |> round(2))
post_df = cbind(post_b_df,post_b0_df,post_b1_df)
rownames(post_df) = rnames
pre_kable = pre_df |>
kable() |>
add_header_above(header = c("Predictor" = 1,
"Odds of Oral Feed \n when oral feed between 0 and 1" = 3,
"Odds of 0% Oral Feed" = 3,
"Odds of 100% Oral Feed" =3)) |>
add_footnote(paste("Posterior variance is estimated to be ", pre_mean[22] |> exp() |> round(2),
" (",pre_lb[22] |> exp() |> round(2),",", pre_ub[22]|> exp()|> round(2),")"))
post_kable = post_df |>
kable() |>
add_header_above(header = c("Predictor" = 1,
"Odds of Oral Feed \n when oral feed between 0 and 1" = 3,
"Odds of 0% Oral Feed" = 3,
"Odds of 100% Oral Feed" =3)) |>
add_footnote(paste("Posterior variance is estimated to be ", post_mean[22] |> exp() |> round(2),
" (",post_lb[22] |> exp() |> round(2),",", post_ub[22]|> exp()|> round(2),")"))Code
pre_kable| Mean | 2.5% | 97.5% | Mean | 2.5% | 97.5% | Mean | 2.5% | 97.5% | |
|---|---|---|---|---|---|---|---|---|---|
| Baseline | 0.96 | 0.21 | 4.47 | 0.20 | 0.02 | 2.32 | 0.02 | 0.00 | 1.11 |
| Attention Score | 1.02 | 0.71 | 1.47 | 0.74 | 0.37 | 1.35 | 1.74 | 0.63 | 5.34 |
| Length of Intubation (d) | 0.86 | 0.76 | 0.98 | 1.35 | 1.13 | 1.67 | 0.69 | 0.45 | 0.99 |
| Single Ventricle w/ Arch Obstruction | 1.57 | 0.71 | 3.36 | 6.19 | 1.90 | 22.14 | 0.76 | 0.03 | 10.62 |
| Two Ventricles w/ Arch Obstruction | 0.81 | 0.42 | 1.59 | 3.43 | 1.11 | 11.58 | 0.52 | 0.08 | 3.05 |
| Age | 0.98 | 0.92 | 1.04 | 0.97 | 0.87 | 1.07 | 1.19 | 1.02 | 1.43 |
| Female | 0.65 | 0.37 | 1.13 | 0.56 | 0.21 | 1.44 | 1.94 | 0.40 | 9.96 |
| a Posterior variance is estimated to be 3 ( 2.08 , 4.22 ) |
Code
post_kable| Mean | 2.5% | 97.5% | Mean | 2.5% | 97.5% | Mean | 2.5% | 97.5% | |
|---|---|---|---|---|---|---|---|---|---|
| Baseline | 0.39 | 0.04 | 3.29 | 0.03 | 0.00 | 0.67 | 0.02 | 0.00 | 37.09 |
| Attention Score | 1.20 | 0.82 | 1.75 | 1.19 | 0.71 | 2.01 | 1.33 | 0.34 | 4.09 |
| Length of Intubation (d) | 0.88 | 0.77 | 1.00 | 1.34 | 1.12 | 1.65 | 0.72 | 0.47 | 1.03 |
| Single Ventricle w/ Arch Obstruction | 1.72 | 0.77 | 3.71 | 6.20 | 1.94 | 21.52 | 0.84 | 0.02 | 13.07 |
| Two Ventricles w/ Arch Obstruction | 0.85 | 0.44 | 1.64 | 3.16 | 1.06 | 9.91 | 0.63 | 0.09 | 3.92 |
| Age | 0.98 | 0.92 | 1.05 | 0.97 | 0.87 | 1.07 | 1.21 | 1.02 | 1.45 |
| Female | 0.64 | 0.36 | 1.11 | 0.54 | 0.20 | 1.40 | 2.47 | 0.46 | 14.76 |
| a Posterior variance is estimated to be 3.07 ( 2.13 , 4.34 ) |
Code
pre_op_pred = readRDS("data/pre_op_pred.rds")
pooled_pre_op_pred =pre_op_pred |> runjags::combine.mcmc(collapse.chains = TRUE)
set.seed(1262023)
sample_pred =sample_n(pooled_pre_op_pred |> as.data.frame(),9)
par(mfrow=c(3,3))
for(i in 1:9) sample_pred[i,] |> as.numeric() |>hist(main = "",
xlab = "")
mtext("Posterior Predictive Samples, Pre-Op Model", side = 3, line = -2, outer = TRUE)
dat[[1]]$Percent_of_feeds_taken_by_mouth_at_discharge |>
hist(main = "Histogram of Percent Oral Feeds",xlab = "Percent Oral Feeds")
Code
# CODE HISTOGRAM
post_op_pred = readRDS("data/post_op_pred.rds")
pooled_post_op_pred =post_op_pred |> runjags::combine.mcmc(collapse.chains = TRUE)
set.seed(1262023)
sample_pred =sample_n(pooled_post_op_pred|> as.data.frame(),9)
par(mfrow=c(3,3))
for(i in 1:9) sample_pred[i,] |> as.numeric() |>hist(main = "",
xlab = "")
mtext("Posterior Predictive Samples, Post-op Model", side = 3, line = -2, outer = TRUE)
### GET ZOIB MODEL FOR DICHOTOMOUS MODEL
# Pre-operation
set.seed(11282023)
tictoc::tic()
# Define a function for model fitting
fit_model <- function(d) {
zoib(
Percent_of_feeds_taken_by_mouth_at_discharge ~
Pre_op_attention_2cat +
Length_of_intubation_days +
Cardiac_Anatomy +
Age_at_Surgery_days +
Female #x1 design matrix
| 1 | #x2 design matrix
Pre_op_attention_2cat +
Length_of_intubation_days +
Cardiac_Anatomy +
Age_at_Surgery_days +
Female | #X3 design matrix
Pre_op_attention_2cat +
Length_of_intubation_days +
Cardiac_Anatomy +
Age_at_Surgery_days +
Female, #x4 design matrix
data = d,
n.response = 1,
zero.inflation = TRUE,
one.inflation = TRUE,
link.mu = "logit",
link.x0 = "logit",
link.x1 = "logit",
random = 0,
n.chain = 4,
n.iter = 3000,
n.thin = 2,
n.burn = 200,
seeds = c(11, 29, 20, 23)
)
}
pre_model_results <- fit_model(dat)
# Save results
saveRDS(pre_model_results, "cat_pre_op_model.rds")
tictoc::toc()
# Post- Operation
set.seed(11282023)
tictoc::tic()
# Define a function for model fitting
fit_model <- function(d) {
zoib(
Percent_of_feeds_taken_by_mouth_at_discharge ~
Post_op_attention_2cat +
Length_of_intubation_days +
Cardiac_Anatomy +
Age_at_Surgery_days +
Female #x1 design matrix
| 1 | #x2 design matrix
Post_op_attention_2cat +
Length_of_intubation_days +
Cardiac_Anatomy +
Age_at_Surgery_days +
Female | #X3 design matrix
Post_op_attention_2cat +
Length_of_intubation_days +
Cardiac_Anatomy +
Age_at_Surgery_days +
Female, #x4 design matrix
data = d,
n.response = 1,
zero.inflation = TRUE,
one.inflation = TRUE,
link.mu = "logit",
link.x0 = "logit",
link.x1 = "logit",
random = 0,
n.chain = 4,
n.iter = 3000,
n.thin = 2,
n.burn = 200,
seeds = c(11, 29, 20, 23)
)
}
# Parallelize model fitting
post_model_results <- dat |> fit_model()
# Save results
saveRDS(post_model_results, "cat_post_op_model.rds")
tictoc::toc()Code
post_model_results = readRDS("data/cat_post_op_model.rds")
pre_model_results = readRDS("data/cat_pre_op_model.rds")
summary_pre_op = pre_model_results$coeff |> summary()
summary_post_op = post_model_results$coeff |> summary()
rnames = c("Baseline","Attention Score","Length of Intubation (d)",
"Single Ventricle w/ Arch Obstruction","Two Ventricles w/ Arch Obstruction",
"Age", "Female")
pre_mean = summary_pre_op$statistics[,"Mean"]
pre_lb = summary_pre_op$quantiles[,"2.5%"]
pre_ub = summary_pre_op$quantiles[,"97.5%"]
post_mean = summary_post_op$statistics[,"Mean"]
post_lb = summary_post_op$quantiles[,"2.5%"]
post_ub = summary_post_op$quantiles[,"97.5%"]
pre_b_df= cbind("Mean"= pre_mean[1:7] |> exp() |> round(2),
"2.5%"= pre_lb[1:7] |> exp() |> round(2),
"97.5%" = pre_ub[1:7] |> exp() |> round(2))
pre_b0_df= cbind("Mean"= pre_mean[8:14] |> exp() |> round(2),
"2.5%"= pre_lb[8:14] |> exp() |> round(2),
"97.5%" = pre_ub[8:14] |> exp() |> round(2))
pre_b1_df= cbind("Mean"= pre_mean[15:21] |> exp() |> round(2),
"2.5%"= pre_lb[15:21] |> exp() |> round(2),
"97.5%" = pre_ub[15:21] |> exp() |> round(2))
pre_df = cbind(pre_b_df,pre_b0_df,pre_b1_df)
rownames(pre_df) = rnames
post_b_df= cbind("Mean"= post_mean[1:7] |> exp() |> round(2),
"2.5%"= post_lb[1:7] |> exp() |> round(2),
"97.5%" = post_ub[1:7] |> exp() |> round(2))
post_b0_df= cbind("Mean"= post_mean[8:14] |> exp() |> round(2),
"2.5%"= post_lb[8:14] |> exp() |> round(2),
"97.5%" = post_ub[8:14] |> exp() |> round(2))
post_b1_df= cbind("Mean "= post_mean[15:21] |> exp() |> round(2),
"2.5%"= post_lb[15:21] |> exp() |> round(2),
"97.5%" = post_ub[15:21] |> exp() |> round(2))
post_df = cbind(post_b_df,post_b0_df,post_b1_df)
rownames(post_df) = rnames
pre_kable = pre_df |>
kable() |>
add_header_above(header = c("Predictor" = 1,
"Odds of Oral Feed \n when oral feed between 0 and 1" = 3,
"Odds of 0% Oral Feed" = 3,
"Odds of 100% Oral Feed" =3)) |>
add_footnote(paste("Posterior variance is estimated to be ", pre_mean[22] |> exp() |> round(2),
" (",pre_lb[22] |> exp() |> round(2),",", pre_ub[22]|> exp()|> round(2),")"))
post_kable = post_df |>
kable() |>
add_header_above(header = c("Predictor" = 1,
"Odds of Oral Feed \n when oral feed between 0 and 1" = 3,
"Odds of 0% Oral Feed" = 3,
"Odds of 100% Oral Feed" =3)) |>
add_footnote(paste("Posterior variance is estimated to be ", post_mean[22] |> exp() |> round(2),
" (",post_lb[22] |> exp() |> round(2),",", post_ub[22]|> exp()|> round(2),")"))Code
pre_kable| Mean | 2.5% | 97.5% | Mean | 2.5% | 97.5% | Mean | 2.5% | 97.5% | |
|---|---|---|---|---|---|---|---|---|---|
| Baseline | 0.86 | 0.35 | 2.22 | 0.08 | 0.01 | 0.38 | 0.05 | 0.00 | 0.65 |
| Attention Score | 1.68 | 0.95 | 2.93 | 0.79 | 0.29 | 2.02 | 5.15 | 1.02 | 32.39 |
| Length of Intubation (d) | 0.87 | 0.76 | 0.98 | 1.35 | 1.13 | 1.67 | 0.72 | 0.49 | 1.03 |
| Single Ventricle w/ Arch Obstruction | 1.59 | 0.74 | 3.32 | 6.54 | 2.08 | 22.83 | 0.66 | 0.02 | 9.69 |
| Two Ventricles w/ Arch Obstruction | 0.85 | 0.44 | 1.62 | 3.23 | 1.07 | 10.50 | 0.66 | 0.10 | 4.08 |
| Age | 0.97 | 0.91 | 1.03 | 0.97 | 0.87 | 1.07 | 1.18 | 1.01 | 1.40 |
| Female | 0.65 | 0.38 | 1.12 | 0.50 | 0.18 | 1.29 | 1.87 | 0.36 | 10.29 |
| a Posterior variance is estimated to be 3.18 ( 2.22 , 4.49 ) |
Code
post_kable| Mean | 2.5% | 97.5% | Mean | 2.5% | 97.5% | Mean | 2.5% | 97.5% | |
|---|---|---|---|---|---|---|---|---|---|
| Baseline | 0.81 | 0.30 | 2.19 | 0.07 | 0.01 | 0.34 | 0.13 | 0.01 | 1.64 |
| Attention Score | 1.45 | 0.79 | 2.54 | 1.09 | 0.44 | 2.67 | 0.72 | 0.11 | 3.91 |
| Length of Intubation (d) | 0.86 | 0.75 | 0.98 | 1.36 | 1.13 | 1.68 | 0.70 | 0.46 | 0.99 |
| Single Ventricle w/ Arch Obstruction | 1.82 | 0.82 | 4.01 | 6.38 | 2.00 | 21.45 | 0.75 | 0.03 | 10.18 |
| Two Ventricles w/ Arch Obstruction | 0.84 | 0.43 | 1.65 | 3.22 | 1.08 | 10.42 | 0.59 | 0.09 | 3.73 |
| Age | 0.99 | 0.93 | 1.05 | 0.97 | 0.87 | 1.07 | 1.19 | 1.02 | 1.42 |
| Female | 0.64 | 0.37 | 1.11 | 0.50 | 0.18 | 1.26 | 1.97 | 0.41 | 9.71 |
| a Posterior variance is estimated to be 3.04 ( 2.1 , 4.27 ) |
3 Appendix
3.1 Check Convergence of MCMC Chains
3.1.1 Traceplots
3.1.1.1 Pre-op 20 imputations model
Code
pooled_pre_op |> traceplot()
3.1.1.2 Post-op 20 imputations model
Code
pooled_post_op |> traceplot()
3.1.2 Autocorrelation plots
3.1.2.1 Pre-op 20 imputations model
Code
pooled_pre_op |> autocorr.plot()
3.1.2.2 Post-op 20 imputations model
Code
pooled_post_op |> autocorr.plot()
3.1.3 Linearity Check
Code
post_pred <- readRDS("data/post_op_pred.rds")
pre_pred <- readRDS("data/pre_op_pred.rds")
post_res <- readRDS("data/post_op_res.rds")
pre_res <- readRDS("data/pre_op_res.rds")
post_pdf <- rep(0,115)
post_rdf <- rep(0,115)
for (i in 1:20){
preds <- rep(0,115)
resids <- rep(0,115)
prd_df1 <- cbind(post_pred[[i]][[1]])
prd_df2 <- cbind(post_pred[[i]][[2]])
prd_df3 <- cbind(post_pred[[i]][[3]])
prd_df4 <- cbind(post_pred[[i]][[4]])
res_df1 <- cbind(post_res[[i]][[1]])
res_df2 <- cbind(post_res[[i]][[2]])
res_df3 <- cbind(post_res[[i]][[3]])
res_df4 <- cbind(post_res[[i]][[4]])
for (j in 1:115){
prd_col1 <- prd_df1[,j]
prd_col2 <- prd_df2[,j]
prd_col3 <- prd_df3[,j]
prd_col4 <- prd_df4[,j]
res_col1 <- res_df1[,j]
res_col2 <- res_df2[,j]
res_col3 <- res_df3[,j]
res_col4 <- res_df4[,j]
preds[j] <- mean(c(prd_col1,prd_col2,prd_col3,prd_col4))
#print(i,j,preds[j])
resids[j] <- mean(c(res_col1,res_col2,res_col3,res_col4))
}
post_pdf <- cbind(post_pdf,preds)
post_rdf <- cbind(post_rdf,resids)
}
# Post operation attention model linearity plots
for (i in 1:20){
pst_df <- cbind(pred=post_pdf[,i+1],resid=post_rdf[,i+1])
pst_df <- as.data.frame(pst_df)
post_plot <- ggplot(pst_df,aes(x=resid,y=pred))+
geom_point(alpha=0.7) +
labs(x="Residuals",y='Predicted Outcome')
disp <- post_plot +
stat_smooth(geom='smooth',method='lm',formula=y~splines::ns(x,knots=c())) +
ggtitle(glue("Linearity of Residuals Post Attention Model Imputed Dataset {i}"))+
theme(plot.title = element_text(face='bold',hjust = 0.5))
print(disp)
}
Code
pre_pdf <- rep(0,115)
pre_rdf <- rep(0,115)
for (i in 1:20){
preds <- rep(0,115)
resids <- rep(0,115)
prd_df1 <- cbind(pre_pred[[i]][[1]])
prd_df2 <- cbind(pre_pred[[i]][[2]])
prd_df3 <- cbind(pre_pred[[i]][[3]])
prd_df4 <- cbind(pre_pred[[i]][[4]])
res_df1 <- cbind(pre_res[[i]][[1]])
res_df2 <- cbind(pre_res[[i]][[2]])
res_df3 <- cbind(pre_res[[i]][[3]])
res_df4 <- cbind(pre_res[[i]][[4]])
for (j in 1:115){
prd_col1 <- prd_df1[,j]
prd_col2 <- prd_df2[,j]
prd_col3 <- prd_df3[,j]
prd_col4 <- prd_df4[,j]
res_col1 <- res_df1[,j]
res_col2 <- res_df2[,j]
res_col3 <- res_df3[,j]
res_col4 <- res_df4[,j]
preds[j] <- mean(c(prd_col1,prd_col2,prd_col3,prd_col4))
#print(i,j,preds[j])
resids[j] <- mean(c(res_col1,res_col2,res_col3,res_col4))
}
pre_pdf <- cbind(pre_pdf,preds)
pre_rdf <- cbind(pre_rdf,resids)
}
# Pre-operation attention model linearity plots
for (i in 1:20){
pst_df <- cbind(pred=post_pdf[,i+1],resid=post_rdf[,i+1])
pst_df <- as.data.frame(pst_df)
post_plot <- ggplot(pst_df,aes(x=resid,y=pred))+
geom_point(alpha=0.7) +
labs(x="Residuals",y='Predicted Outcome')
disp <- post_plot +
stat_smooth(geom='smooth',method='lm',formula=y~splines::ns(x,knots=c())) +
ggtitle(glue("Linearity of Residuals Pre Attention Model Imputed Dataset {i}"))+
theme(plot.title = element_text(face='bold',hjust = 0.5))
print(disp)
}
3.1.4 Posterior Predictive Check
Code
pre_op_pred = readRDS("data/pre_op_pred.rds")
pooled_pre_op_pred =pre_op_pred |> runjags::combine.mcmc(collapse.chains = TRUE)
post_op_pred = readRDS("data/post_op_pred.rds")
pooled_post_op_pred =post_op_pred |> runjags::combine.mcmc(collapse.chains = TRUE)
#| cache: true
pre_sample_pred =sample_n(pooled_pre_op_pred|> as.data.frame(),1000)
post_sample_pred = sample_n(pooled_post_op_pred|> as.data.frame(),1000)
pre_q25 = rep(NA, 1000)
pre_q975 = rep(NA, 1000)
post_q25 = rep(NA, 1000)
post_q975 = rep(NA,1000)
for(i in 1:1000){
m = pre_sample_pred[i,] |> as.numeric()
pre_q25[i] = quantile(m, .025)
pre_q975[i] = quantile(m, .975)
m = post_sample_pred[i,] |> as.numeric()
post_q25[i] = quantile(m, .025)
post_q975[i] = quantile(m, .975)
}
q25 = dat[[1]]$Percent_of_feeds_taken_by_mouth_at_discharge |>
quantile(.025)
q975 = dat[[1]]$Percent_of_feeds_taken_by_mouth_at_discharge |>
quantile(.975)
post_q25 |> hist(main = "Histogram of 2.5th Percentile of Regenerated Samples \n Post-op Model")
abline(v = q25, col = "blue")
Code
pre_q25 |> hist(main = "Histogram of 2.5th Percentile of Regenerated Samples \n Pre-op Model")
abline(v = q25, col = "blue")
Code
post_q975 |> hist(main = "Histogram of 97.5th Percentile of Regenerated Samples \n Post-op Model")
abline(v = q975, col = "blue")
Code
pre_q975 |> hist(main = "Histogram of 97.5th Percentile of Regenerated Samples \n Pre-op Model")
abline(v = q975, col = "blue")
Code
d = dat[[1]]$Percent_of_feeds_taken_by_mouth_at_discharge
d_no01 = d[d !=0 & d != 1]
d_no01 |> hist(main = "Histogram of Percentage Oral Feed \n excluding 0 and 1")
Code
post_q25[d !=0 & d != 1] |> hist(main = "Histogram of 2.5th Percentile of Regenerated Samples \n Post-op Model")
abline(v = quantile(d_no01,.025), col = "blue")
Code
pre_q25[d !=0 & d != 1] |> hist(main = "Histogram of 2.5th Percentile of Regenerated Samples \n Pre-op Model")
abline(v = quantile(d_no01,.025), col = "blue")
Code
post_q975[d !=0 & d != 1] |> hist(main = "Histogram of 97.5th Percentile of Regenerated Samples \n Post-op Model")
abline(v = quantile(d_no01,.975), col = "blue")
Code
pre_q975[d !=0 & d != 1] |> hist(main = "Histogram of 97.5th Percentile of Regenerated Samples \n Pre-op Model")
abline(v = quantile(d_no01,.975), col = "blue")
Code
d = dat[[1]]$Percent_of_feeds_taken_by_mouth_at_discharge